################################################################################################################## # Functions that implement Individual testing, Master pool testing, Dorfman testing, Array testing ############################################################################# # R function: This function simulates individual level testing and stores the # testing responses in accordance to the data structure # required to fit the Bayesian GT regression model presented # in McMahan et al. via the R function Bayes.reg.GT # # Input: # Y.true = The true statuses of the individuals # Se = The testing sensitivity used for both pools and individual testing # Sp = The testing specificity used for both pools and individual testing # Ind.decode<-function(Y.true,Se,Sp){ N<-length(Y.true) Y<-rbinom(N,1,(Se*Y.true+(1-Sp)*(1-Y.true))) return(Y) } ##################################################################### # R function: This function simulates Initial pool testing and stores the # testing responses in accordance to the data structure # required to fit the Bayesian GT regression model presented # in McMahan et al. via the R function Bayes.reg.GT # # Input: # Y.true = The true statuses of the individuals # Se = A vector of testing sensitivities, where the first element is the # testing sensitivity for the pools and the second entry is the # test sensitivity for individual testing # Sp = A vector of testing specificities, where the first element is the # testing specificity for the pools and the second entry is the # test specificity for individual testing # cj = pool size to be used (Note: The number of individuals N should be # evenly divisible by cj, this is only for decoding purposes; i.e., # the regression methods do not require this condition) Pool.test<-function(Y.true,Se,Sp,cj){ N<-length(Y.true) Jmax<-N/cj J<-1 Y<-matrix(-99,nrow=N, ncol=3) Z<-matrix(-99,nrow=Jmax,ncol=cj+3) for(j in 1:(N/cj)){ prob<-ifelse(sum(Y.true[((j-1)*cj+1):(j*cj)])>0,Se,1-Sp) Z[J,1]<-rbinom(n=1,size=1,prob=prob) Z[J,2]<-cj Z[J,3]<-1 Z[J,4:(cj+3)]<-((j-1)*cj+1):(j*cj) Y[((j-1)*cj+1):(j*cj),2]<-1 Y[((j-1)*cj+1):(j*cj),3]<-J J<-J+1 } J<-J-1 Z<-Z[1:J,] return(list("Z"=Z, "Y"=Y)) } ##################################################################### # R function: This function simulates Dorfman decoding and stores the # testing responses in accordance to the data structure # required to fit the Bayesian GT regression model presented # in McMahan et al. via the R function Bayes.reg.GT # # Input: # Y.true = The true statuses of the individuals # Se = A vector of testing sensitivities, where the first element is the # testing sensitivity for the pools and the second entry is the # test sensitivity for individual testing # Sp = A vector of testing specificities, where the first element is the # testing specificity for the pools and the second entry is the # test specificity for individual testing # cj = pool size to be used (Note: The number of individuals N should be # evenly divisible by cj, this is only for decoding purposes; i.e., # the regression methods do not require this condition) Dorfman.decode.diff.error<-function(Y.true,Se,Sp,cj){ N<-length(Y.true) Jmax<-N+N/cj J<-1 Y<-matrix(-99,nrow=N, ncol=4) Z<-matrix(-99,nrow=Jmax,ncol=cj+3) for(j in 1:(N/cj)){ prob<-ifelse(sum(Y.true[((j-1)*cj+1):(j*cj)])>0,Se[1],1-Sp[1]) Z[J,1]<-rbinom(n=1,size=1,prob=prob) Z[J,2]<-cj Z[J,3]<-1 Z[J,4:(cj+3)]<-((j-1)*cj+1):(j*cj) Y[((j-1)*cj+1):(j*cj),2]<-1 Y[((j-1)*cj+1):(j*cj),3]<-J J<-J+1 if(Z[J-1,1]==1){ for(k in ((j-1)*cj+1):(j*cj)){ prob<-ifelse(Y.true[k]>0,Se[2],1-Sp[2]) Z[J,1]<- rbinom(n=1,size=1,prob=prob) Z[J,2]<-1 Z[J,3]<-2 Z[J,4]<-k Y[k,2]<-2 Y[k,4]<-J J<-J+1 } } } J<-J-1 Z<-Z[1:J,] return(list("Z"=Z, "Y"=Y)) } ##################################################################### # R function: This function simulates Dorfman decoding and stores the # testing responses in accordance to the data structure # required to fit the Bayesian GT regression model presented # in McMahan et al. via the R function Bayes.reg.GT # # Input: # Y.true = The true statuses of the individuals # Se = The testing sensitivity used for both pools and individual testing # Sp = The testing specificity used for both pools and individual testing # cj = pool size to be used (Note: The number of individuals N should be # evenly divisible by cj, this is only for decoding purposes; i.e., # the regression methods do not require this condition) Dorfman.decode.same.error<-function(Y.true,Se,Sp,cj){ N<-length(Y.true) Jmax<-N+N/cj J<-1 Y<-matrix(-99,nrow=N, ncol=4) Z<-matrix(-99,nrow=Jmax,ncol=cj+3) for(j in 1:(N/cj)){ prob<-ifelse(sum(Y.true[((j-1)*cj+1):(j*cj)])>0,Se,1-Sp) Z[J,1]<-rbinom(n=1,size=1,prob=prob) Z[J,2]<-cj Z[J,3]<-1 Z[J,4:(cj+3)]<-((j-1)*cj+1):(j*cj) Y[((j-1)*cj+1):(j*cj),2]<-1 Y[((j-1)*cj+1):(j*cj),3]<-J J<-J+1 if(Z[J-1,1]==1){ for(k in ((j-1)*cj+1):(j*cj)){ prob<-ifelse(Y.true[k]>0,Se,1-Sp) Z[J,1]<- rbinom(n=1,size=1,prob=prob) Z[J,2]<-1 Z[J,3]<-1 Z[J,4]<-k Y[k,2]<-2 Y[k,4]<-J J<-J+1 } } } J<-J-1 Z<-Z[1:J,] return(list("Z"=Z, "Y"=Y)) } ##################################################################### # R function: This function simulates Array decoding and stores the # testing responses in accordance to the data structure # required to fit the Bayesian GT regression model presented # in McMahan et al. via the R function Bayes.reg.GT # # Input: # Y.true = The true statuses of the individuals # Se = A vector of testing sensitivities, where the first element is the # testing sensitivity for the pools and the second entry is the # test sensitivity for individual testing # Sp = A vector of testing specificities, where the first element is the # testing specificity for the pools and the second entry is the # test specificity for individual testing # cj = row and column pool sizes to be used (Note: The number of individuals # N should be evenly divisible by cj^2, this is only for decoding # purposes; i.e., the regression methods do not require this condition) Array.decode.diff.error<-function(Y.true, Se, Sp, cj){ N<-length(Y.true) Jmax<-2*N/cj +N J<-1 AT<-N/(cj^2) Y<-matrix(-99,nrow=N, ncol=5) Z<-matrix(-99,nrow=Jmax,ncol=cj+3) Y.A<-array(-99,c(cj,cj,AT)) ID.A<-array(-99,c(cj,cj,AT)) ind<-1 for(i in 1:AT){ for(j in 1:cj){ for(m in 1:cj){ Y.A[m,j,i]<-Y.true[ind] ID.A[m,j,i]<-ind ind<-ind+1 }}} for(s in 1:AT){ array.yk<-Y.A[,,s] array.id<-ID.A[,,s] a<-rep(0,nrow(array.yk)) b<-rep(0,ncol(array.yk)) for(i in 1:cj){ prob<- ifelse(sum(array.yk[i,])>0, Se[1], 1-Sp[1]) g<- rbinom(n=1,size=1,prob=prob) a[i]<-g Z[J,1]<-g Z[J,2]<-cj Z[J,3]<-1 Z[J,4:(cj+3)]<-array.id[i,] Y[array.id[i,],2]<-2 Y[array.id[i,],3]<-J J<-J+1 } for(j in 1:cj){ prob<- ifelse(sum(array.yk[,j])>0, Se[1], 1-Sp[1]) g<- rbinom(n=1,size=1,prob=prob) b[j]<-g Z[J,1]<-g Z[J,2]<-cj Z[J,3]<-1 Z[J,4:(cj+3)]<-array.id[,j] Y[array.id[,j],4]<-J J<-J+1 } if(sum(a)>0 & sum(b)>0){ array.yk1<-as.matrix(array.yk[(a==1),(b==1)]) array.id1<-as.matrix(array.id[(a==1),(b==1)]) for(i in 1:nrow(array.yk1)){ for(j in 1:ncol(array.yk1)){ prob<- ifelse(array.yk1[i,j]>0, Se[2], 1-Sp[2]) g<- rbinom(n=1,size=1,prob=prob) Z[J,1]<-g Z[J,2]<-1 Z[J,3]<-2 Z[J,4]<-array.id1[i,j] Y[array.id1[i,j],2]<-3 Y[array.id1[i,j],5]<-J J<-J+1 }}} if(sum(a)>0 & sum(b)==0){ array.yk1<-as.matrix(array.yk[(a==1),]) array.id1<-as.matrix(array.id[(a==1),]) for(i in 1:nrow(array.yk1)){ for(j in 1:ncol(array.yk1)){ prob<- ifelse(array.yk1[i,j]>0, Se[2], 1-Sp[2]) g<- rbinom(n=1,size=1,prob=prob) Z[J,1]<-g Z[J,2]<-1 Z[J,3]<-2 Z[J,4]<-array.id1[i,j] Y[array.id1[i,j],2]<-3 Y[array.id1[i,j],5]<-J J<-J+1 }}} if(sum(a)==0 & sum(b)>0){ array.yk1<-as.matrix(array.yk[,(b==1)]) array.id1<-as.matrix(array.id[,(b==1)]) for(i in 1:nrow(array.yk1)){ for(j in 1:ncol(array.yk1)){ prob<- ifelse(array.yk1[i,j]>0, Se[2], 1-Sp[2]) g<- rbinom(n=1,size=1,prob=prob) Z[J,1]<-g Z[J,2]<-1 Z[J,3]<-2 Z[J,4]<-array.id1[i,j] Y[array.id1[i,j],2]<-3 Y[array.id1[i,j],5]<-J J<-J+1 }}} } J<-J-1 Z<-Z[1:J,] return(list("Z"=Z, "Y"=Y)) } #End function ##################################################################### # R function: This function simulates Array decoding and stores the # testing responses in accordance to the data structure # required to fit the Bayesian GT regression model presented # in McMahan et al. via the R function Bayes.reg.GT # # Input: # Y.true = The true statuses of the individuals # Se = The testing sensitivity used for both pools and individual testing # Sp = The testing specificity used for both pools and individual testing # cj = row and column pool sizes to be used (Note: The number of individuals # N should be evenly divisible by cj^2, this is only for decoding # purposes; i.e., the regression methods do not require this condition) Array.decode.same.error<-function(Y.true, Se, Sp, cj){ N<-length(Y.true) Jmax<-2*N/cj +N J<-1 AT<-N/(cj^2) Y<-matrix(-99,nrow=N, ncol=5) Z<-matrix(-99,nrow=Jmax,ncol=cj+3) Y.A<-array(-99,c(cj,cj,AT)) ID.A<-array(-99,c(cj,cj,AT)) ind<-1 for(i in 1:AT){ for(j in 1:cj){ for(m in 1:cj){ Y.A[m,j,i]<-Y.true[ind] ID.A[m,j,i]<-ind ind<-ind+1 }}} for(s in 1:AT){ array.yk<-Y.A[,,s] array.id<-ID.A[,,s] a<-rep(0,nrow(array.yk)) b<-rep(0,ncol(array.yk)) for(i in 1:cj){ prob<- ifelse(sum(array.yk[i,])>0, Se, 1-Sp) g<- rbinom(n=1,size=1,prob=prob) a[i]<-g Z[J,1]<-g Z[J,2]<-cj Z[J,3]<-1 Z[J,4:(cj+3)]<-array.id[i,] Y[array.id[i,],2]<-2 Y[array.id[i,],3]<-J J<-J+1 } for(j in 1:cj){ prob<- ifelse(sum(array.yk[,j])>0, Se, 1-Sp) g<- rbinom(n=1,size=1,prob=prob) b[j]<-g Z[J,1]<-g Z[J,2]<-cj Z[J,3]<-1 Z[J,4:(cj+3)]<-array.id[,j] Y[array.id[,j],4]<-J J<-J+1 } if(sum(a)>0 & sum(b)>0){ array.yk1<-as.matrix(array.yk[(a==1),(b==1)]) array.id1<-as.matrix(array.id[(a==1),(b==1)]) for(i in 1:nrow(array.yk1)){ for(j in 1:ncol(array.yk1)){ prob<- ifelse(array.yk1[i,j]>0, Se, 1-Sp) g<- rbinom(n=1,size=1,prob=prob) Z[J,1]<-g Z[J,2]<-1 Z[J,3]<-1 Z[J,4]<-array.id1[i,j] Y[array.id1[i,j],2]<-3 Y[array.id1[i,j],5]<-J J<-J+1 }}} if(sum(a)>0 & sum(b)==0){ array.yk1<-as.matrix(array.yk[(a==1),]) array.id1<-as.matrix(array.id[(a==1),]) for(i in 1:nrow(array.yk1)){ for(j in 1:ncol(array.yk1)){ prob<- ifelse(array.yk1[i,j]>0, Se, 1-Sp) g<- rbinom(n=1,size=1,prob=prob) Z[J,1]<-g Z[J,2]<-1 Z[J,3]<-1 Z[J,4]<-array.id1[i,j] Y[array.id1[i,j],2]<-3 Y[array.id1[i,j],5]<-J J<-J+1 }}} if(sum(a)==0 & sum(b)>0){ array.yk1<-as.matrix(array.yk[,(b==1)]) array.id1<-as.matrix(array.id[,(b==1)]) for(i in 1:nrow(array.yk1)){ for(j in 1:ncol(array.yk1)){ prob<- ifelse(array.yk1[i,j]>0, Se, 1-Sp) g<- rbinom(n=1,size=1,prob=prob) Z[J,1]<-g Z[J,2]<-1 Z[J,3]<-1 Z[J,4]<-array.id1[i,j] Y[array.id1[i,j],2]<-3 Y[array.id1[i,j],5]<-J J<-J+1 }}} } J<-J-1 Z<-Z[1:J,] return(list("Z"=Z, "Y"=Y)) } #End function ########################################### # Function used to compute the coverage # probabilities cov.func<-function(data,sim,P,alpha,true){ cov<-matrix(-99,sim,P) for(s in 1:sims){ for(p in 1:P){ int<-quantile(data[,p,s],c(alpha/2,1-alpha/2)) cov[s,p]<-int[1]<=true[p] & true[p]<=int[2] }} cov<-apply(cov,2,mean) return(cov) }